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O ' Abstract 

From a variety of long simulations of all-atom butane using both stochastic and 

fully-solved molecular dynamics, we have uncovered striking generic behavior which 

also occurs in one-dimensional systems. We find an apparently universal distribution 

of transition event durations, as well as a characteristic speed profile along the reaction 

coordinate. An approximate analytic distribution of event durations, derived from 

a one-dimensional model, correctly predicts the asymptotic behavior of the universal 

£> ■ distribution for both short and long durations. 
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1 Introduction 

Ensuring appropriate behavior in molecular simulations has long been a topic of interest, 
and butane has been an important test molecule as the simplest alkane with a dihedral 
degree of freedom. In early work on united-atom models Rebertus, Berne and Chandler 
considered solvent effects on butane conformational equilibrium |I| , while Levy, Karplus and 
McCammon studied stochastic butane dynamics 0. Pastor and Karplus later investigated 
the importance of inertial effects in Langevin dynamics, also using a united atom model of 
butane || . Tobias and Brooks discussed the importance of including hydrogen atoms when 
considering solvent effects ||. 

One focus of molecular studies has been the dependence of transition rates on the choice 
of dynamics. In a study of ethylene glycol, which is polar but also possesses a single dihedral 
angle, Widmalm and Pastor compared rates and other timescales in molecular and Langevin 
dynamics ||. Similarly, Loncharich, Brooks and Pastor explored the effects of the Langevin 
friction constant on transitions rates in the larger "NAMA" molecule with an eye toward 
speeding up the generation of conformational ensembles in peptides and proteins |J . 

Work in one-dimensional systems, in addition to considering the effects of stochastic 
dynamics on reaction rates (e.g., |7|, ||, |9], [U], [TTJ, |12|1), has explored transition phenomena on 



still-shorter "microscopic" timescales. Dykman and coworkers have extensively investigated 
the effects of different noise types and non-equilibrium driving forces on transition events 



themselves in stochastic settings |13L 14, [Lq, [Lq, O, M. The present authors previously 



considered the influence of the local, non-uniform speed along the reaction coordinate on a 



"dynamic importance sampling" approach to rate computations [[HJ . 

Many "microscopic" questions have remained unanswered regarding molecular transition 
events. What is their duration? Is the reaction coordinate traversed at a uniform speed, 
or is progress more rapid in some regions? Is there sensitivity to the choice of dynamics 
or potential? The need for detailed answers is more pressing in molecular systems because 



dynamic importance sampling [£(J, [H], [ij|] and "transition path sampling" methods |22], ^3[ 
for computing reaction rates require ensembles of transition event trajectories. A quantitative 
understanding of transition events could accelerate ensemble generation in these approaches, 
making larger systems amenable to study. Moreover, because ensembles of transition events 
are required for rate estimates, the durations set a fundamental minimum on the computer 



time required to perform such calculations [HI . 

The present study provides answers to these questions, using a spectrum of simulations. 
First, using the CHARMM potential, all-atom butane is studied in explicit water and acetic 
acid. Low- and high-friction Langevin simulations are also performed with CHARMM, along 
with a separate high-friction run using the AMBER potential. We compare the butane event- 
duration distributions and reaction-coordinate speed profiles with those for one dimensional 
stochastic systems with white and colored noise. 

Two interesting generic features emerge. First, after a simple linear re-scaling by the av- 
erage, the distribution of transition event durations appears to be universal across dynamics, 
solvent, and potential type — regardless of whether butane or a one-dimensional system is 
considered. The distribution embodies behaviors predicted analytically from a simple one- 
dimensional model. Second, the reaction-coordinate speed profiles of butane transitions also 
display a common, accelerating trend, which matches that of colored-noise, one- dimensional 
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Figure 1: A single transition event in a butane simulation, illustrating the transition event 
duration, t b . The precise definition of t^ is given in Sec. fO 



simulations. 

In the next section, |^, we discuss simulation models and methods. Section ^| presents an 
analytic prediction for the distribution of event durations. Sec. |] contains the simulation 
results, and we give a summary and some conclusions in Section |5|. 



2 Simulation Models and Set-Up 

Transition behavior was studied in a variety of settings — from all-atom, fully solvated molec- 
ular dynamics to one-dimensional Langevin dynamics. This section describes the protocols 
and parameters used. 

2.1 Fully-Solvated All- Atom Models of Butane 

To examine solvent effects on barrier crossing behavior, and to ensure a good standard 
of comparison for stochastic simulations, we performed two explicit solvent simulations of 
butane. In one case the solvent was water and in the other it was acetic acid. In both 
cases the CHARMM program and parameters was used. The butane molecule was lightly 
restrained to the center of the simulation cell using a harmonic, molecule-specific, center-of- 
mass energy term. The simulations were performed with constant pressure and temperature 
and a non-bonded list cut-off of 12 A. The non-bonded list was generated to 15 A and the van 
der Waals force switched to zero over a 12 A distance. Electrostatics were shifted over the 
full 15 A range to cut-off. Dynamics used the leapfrog algorithm (e.g., H|) in CHARMM, 
which for each particle's position x and velocity v is given by 

x(t+At) = x(t) + Atv(t + |Ai) (1) 

v(£+±Ai) = v(t-|At)+Atf(t)/m, (2) 

where f = — V X C/ for potential energy U and m is the particle's mass. 

The acetic acid simulations consisted of 253 explicit molecules of acetic acid and box 
dimensions of roughly 65 A on a side. This created the appropriate density of solvent for a 
temperature of 300 K. The water simulations contained 861 TIP3 water molecules and used 
a box size of roughly 30 A. 

2.2 All-Atom Butane in a Stochastic "Solvent" 

To maintain as much consistency as possible between the explicit solvent and the implicit 
solvent simulations, the stochastic simulations used similar options to the explicit solvent 
simulations. In particular, the non-bonded energies were cut off at 12 A, and a center-of- 
mass restraint was again used on the butane molecule. Temperature was maintained at 300 
K. As would be expected, no periodic images were used and there is no pressure control in 
these stochastic simulations. Dynamics were performed with the leapfrog Langevin dynamics 
integrator in CHARMM, namely f|, 

x(i + A£) = x(t) + At v(t + ±Ai) (3) 

x(i) -x(i-Ai) 1- \l^t 



v(£+±Ai) = v(t-|At) + 



At 1 + i 7 At 



,. . ,f(t) + R(£) 

+(A ' /m) TTW (4) 



v(«) = ,/l + ±7Ae[v(«+iAe)+v(t-±Ai)]/2, (5) 



where R is a white-noise frictional force, every component of which is chosen from a Gaussian 
distribution of zero mean and variance 2m , ykBT '/ At. We used two values of the friction 
parameter (7 = 5,50 ps _1 ) to implicitly model solvent effects; only the carbon atoms were 
treated stochastically, while the hydrogen dynamics were governed by the ordinary (7 = 0) 
leapfrog Verlet algorithm ([[]). 

An additional Langevin simulation was performed using the AMBER potential within 
the Molecular Modeling Tool Kit (MMTK) |2§. For this "overdamped" ("Brownian") sim- 
ulation, dynamics were governed by 



x(t + At) = x(i) + (f/mrf)At + Ax R 



(6) 



The noise term Ax^ was also chosen to be "white" - in this case, selected independently 
at every time step from a Gaussian of zero mean and variance 



a 2 = 2Atk B T/mj . 



(7) 



All atoms were treated stochastically using the convenient, but we note that overdamped 
dynamics only make predictions relative to an arbitrary overall timescale || embodied here 
in the (convenient) choice 7 = 150. 

2.3 One-Dimensional Models 
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Figure 2: The bistable, one-dimensional potential (§) simulated under a variety of conditions 
discussed in the text. It is shown here for a barrier height of E^ = 7kgT and lengthscale 
1 = 1. 



We also studied the simple bistable potential discussed in |L9[ and many other places, 
namely, 

U(x) = E b [{x/lf-l] 2 , (8) 



where Ef, is the barrier height and I the length scale. This potential is shown in Fig. |2|. All 
one-dimensional simulations used overdamped Langevin dynamics @ with — for simplicity 
- mass, friction, and the thermal energy scale set to unity: m = 7 = fcgT = 1. 

Two types of noise were studied. The first was uncorrelated "white" noise, as dis- 
cussed above. We also employed exponentially-correlated noise simulated via an Ornstein- 
Uhlenbeck process ||26|| , in which Axr = y is considered a position variable executing in- 
dependent overdamped dynamics (^) in a harmonic potential U(y) = (m , y/2T)y 2 subject 
to white noise with fluctuations given by at = (2At/r). This formulation yields the auto- 
correlation function (y(0)y(t)) = <r 2 exp (— t/r), so that the magnitude of the correlated 
noise fluctuations is unchanged from the white-noise case. 



3 Theoretical Analysis of a Simple Model 
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Figure 3: A simple one-dimensional potential used to calculate probabilities of event dura- 
tions, t b . 

This section discusses an extremely idealized potential characterized by a transition- 
event-duration distribution with features common among all of the models studied. The 
reader solely interested in the simulation results could safely pass over this discussion and 
refer to the asymptotic results quoted later. 

The one-dimensional potential is shown in Fig. |3|, and the initial state A is defined to be 
x < 2 and state B is x > 3. The motivation for studying such an unphysical model is based 
on the fundamental observation made in the present report (see next section): because the 
distribution of durations appears to be universal, it should be derivable from even a trivial 
model. 

The analysis is based on the simulation-step probabilities governing the overdamped 
dynamics @ with white noise. Explicitly, the potential of Fig. ^ is given by 

U(x) = (x < 2) 

= (E b /L)(x-2) (2<x<3) (9) 

= E b (x > 3) . 

For convenience, we introduce the notation 5x = fAt/rwy = EbAt/rwyL for the inclined 
region 2 < x/L < 3. Because fluctuation increments Axr are chosen from a Gaussian 
distribution, the probability density for choosing an increment Axj = Xj + i — Xj in the 
inclined region is given by 

T(Axj) = ^^ exp [-(Axj - 5xf/2a 2 ] . (10) 
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The probability density of a trajectory of n& = t b /At steps, ( n ' b = (x , Xi, . . . , x nb ), is simply 
the product of single-step densities: 



n h —\ 

q(C) = n n^ 3 ) , 

3=0 

which, because each of the Gaussian form ([T0|), is properly normalized according to 



;n) 



dCQ(C 



(12) 



where d( rib = n™=o ^Axj. 

The relative probability of a given event duration, £& = rib At, is a product of two factors: 
(i) the total number of probable n&-step trajectories — whether crossing or not — multiplied 
by (ii) the fraction of n^-step trajectories which cross the barrier (i.e., reach the barrier top 
at x > 3; see Fig. |). Factor (i) is proportional to a" 6 , because the fluctuation width in 
(|l0|) measures the extent of probable steps. The second factor, the fraction of successful 
trajectories, may be written formally as 



h 



dC nb Q(C b )h A (xo)h B (x 



n b ) i 



(13) 



where the indicator function hy{x) is unity when x is in state Y and zero otherwise. That is, 
the probability of a transition event occurring in rib = h/ At steps is proportional to cr nb f b . 
We estimate the fraction of successful trajectories, /&, by considering only those tra- 
jectories near to the optimal r^-step crossing trajectory (which differs from the unique, 
overall-optimal trajectory -- for the the optimal n& -- discussed in [[H|). The optimal rib- 
step trajectory for the constant-force potential of Fig. |3] consists of uniform steps of length 
L/rib, and we estimate the number of nearby trajectories as those occurring in a "fluctuation 
volume" a Ub about the optimal; that is, we estimate 



f b oc c tb/At exp 



rwy 



t b (- + ^ 
2k bT \tb rwyLj 



V 



(14) 



where Cq is an unknown constant, presumably of order unity. 

Finally, the estimate for the probability distribution of event durations is proportional 
to the product of o Ub and fb, namely, 



p{t b ) oc { Cl a) tb ' M 



exp 



2ksT \tb m^L ) 



(15) 



where c\ is again an unknown constant. Using the identity a n = e nloga , we can rewrite this 
approximation as a two-parameter form for purposes of understanding its behavior, 



p(tb) oc exp 



-ii, ( Y + d 



(16) 



The asymptotic behavior is of particular interest and is given by 



p(t b ^0) ~ exp(-c 2 /t b ), (17) 

p(tb — ► oo) ~ exp (—d 2 t b ) . (18) 



It is worth noting that all derivatives of p vanish as tb — ► 0. 



4 Results 

4.1 Transition- Event Durations 

We here present results regarding the durations of transition events, if,; see Fig. [I]. The 
universal behavior of the distribution of durations is demonstrated, and limiting behavior is 
discussed by comparison to the analysis of the previous section. The average durations are 
also compared to other timescales which distinguish each system, the correlation time and 
inverse reaction rate. 
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Figure 4: Collapse of distributions of transition event durations, £&, after simple linear rescaling. 
Both butane and one-dimensional simulation data are shown, and the probability density function, 
p, is plotted for the durations normalized by the average of each data set. The lower plot, shows 
identical data on a logarithmic scale to show the common behavior in the tails of the distributions. 
The butane data sets represent water- and acetic-acid-solvated molecular dynamics (MD) with the 
CHARMM potential, as well as low- and high-friction Langevin dynamics (LD) with CHARMM, 
plus a high- friction (overdamped) Langevin simulation using AMBER. In one dimension (fD), 
overdamped Langevin dynamics were used, with both white and colored noise. The empirical fit is 
given in di~9| ) and discussed in the text. 
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The definition of a transition duration is necessarily somewhat arbitrary. Here we used 
the time between the last exit from the initial state and the first entrance to the final state; 
in other words, tb is the amount of time spent continuously by a trajectory in the transition 
region. For transitions between the gauche— and trans states in butane, we chose the the 
transition region to be 100 < <fi < 150; considering symmetry, the region 210 < <p < 260 was 
chosen for gauche+ and trans. Following convention, the trans state is centered at <fi = 180. 

The principal results for the distributions of event durations are embodied in Fig. [|. 
The simple linear rescalings, by the system-specific average durations (t b ), cause the data to 
collapse onto a fairly well-determined universal curve. The empirical fit to the data, detailed 
below, confirms the asymptotic behaviors predicted analytically in Eqs. (|P7D and (|i~8|): the 



gap at small tb is clearly consistent with e _1//f6 behavior (but not with a more common form 
like tb q e~ tb with exponent q > 0), and the logarithmic scale of Fig. |](b) reveals the large tb 
behavior to be a simple exponential decay. 

The empirical fit shown in Fig. |] depicts the form 



p(t b ) =M x exp 



b + t b r h 2 



(19) 



where M is a normalization factor and the three parameters used are b = 1.5, c = 1.4, r = 5. 
The form was designed, by trial and error, both to display the correct asymptotic behaviors 
and to capture the inflection point to the right of the peak visible on the logarithmic plot, 
Fig. 1(b). 

4.2 Timescales 

Another question that can be answered based on our data is whether the crossing event 
durations in a particular system are correlated with other fundamental timescales, like the 
correlation time and the (inverse) reaction rates. Table p] apparently answers this question 
in the negative. Note that, in the table, the dihedral correlation time is defined by 

POO 

W = (0(0)0(0))^ / d£(0(O)0(t)>, (20) 



Table 1: Comparison of timescales for butane dynamics with different explicit and implicit 
solvents. The rates, k, are between the trans (t) and gauche (g) states, and the reciprocals 
define waiting times between transition events. For definitions of the dihedral correlation 
times r corr and Td eC ay, see Eq. (^Up and the subsequent text. The average transition-event 
duration is denoted by (tb). All times are in picoseconds. 

System 1/kg-H l/h-*g 7"decay ^corr (h) 



MD /acetic acid 


55.6 


294 


0.40 


0.02 


0.136 


MD /water 


35.1 


115 


0.21 


0.04 


0.145 


LD/( 7 = 5) 


38.5 


154 


0.35 


0.02 


0.146 


LD/(7 = 50) 


32.8 


151 


0.06 


0.076 


0.200 
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where the angular brackets refer here to time averages within the trans state only. Since the 
dihedral-angle correlation function of the integrand exhibits strong oscillations in molecular 
dynamics simulations, Table |1| also includes the decay time of the exponential "envelope," 
Tdecay The data also suggest that the rates are not simply related to correlation times. 

4.3 Reaction-Coordinate Speed Profiles 



The reaction-coordinate speed profile |L9| depicts the local speed of progress along the re- 
action coordinate, as exemplified in Fig. |5|. This sub-picosecond behavior (see Table [I]) is 
uniquely available from simulations, and critically important for dynamic importance sam- 
pling methods for rate calculations: see Ref. [|19| . 



As in the case of the the crossing time distributions, the results tell a simple and consistent 
story. All of the butane simulations show a characteristic "accelerating" speed profile during 
transitions, as indicated by Fig. [|. Absolute speeds do differ, as is evident from the average 
crossing times of Table |]. The consistent, accelerating profiles from the butane simulations 
contrast sharply with the white-noise profile of a one-dimensional simulation discussed in 
detail in Ref. |19| and also shown in Fig. |6|. 

Because even the Langevin butane simulations — which were performed with white noise 
- exhibit the accelerating behavior, one can also conclude that molecular fluctuations are 
inherently colored (i.e., correlated). The reason is straightforward. If one focuses on the 
dynamics of particular butane atom, say a methyl carbon, then according to Eq. ([]) that 
atom moves based on the force it feels and a white noise increment. The force itself, however, 
reflects the relatively long-timescale motions executed by the whole molecule. Hence the 
molecule acts as its own memory, and correlations are not surprising regardless of the type 
of noise modeled in Axr. 
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Figure 5: The "average crossing speed" during gauche— to trans transitions in various 
model butane systems. We plot the speed (A(j)/At) during transition events only, normal- 
ized by the overall average during these events. The five butane data sets represent: (i) 
water-solvated molecular dynamics with the CHARMM potential; (ii) acetic-acid-solvated 
molecular dynamics with CHARMM; (iii) low-friction Langevin dynamics with CHARMM; 
(iv) high-friction Langevin dynamics with CHARMM; and, (v) high-friction Langevin dy- 
namics with the AMBER potential. The absolute crossing speeds can be inferred from the 
data of Table |l[ For reference, the bottom panel shows the free energy as a function of 
dihedral angle in the transition region. 
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Figure 6: As in Fig. [5], the reaction coordinate speed profile, now for the one- dimensional 
bistable well depicted in Fig. |2|. The average reaction coordinate speeds (Ax/ At) are plotted 
as a function of x for high-friction Langevin simulations, in arbitrary units; the computations, 
including colored and white noise, are described in Sec. p75| . The "accelerating" trend of the 
colored noise data matches the butane speed profiles, while the symmetric white-noise data 
do not. 
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5 Summary and Conclusions 

By comparing butane simulations with explicit and implicit solvents, along with simple 
white- and colored- noise one-dimensional (ID) systems, we have found a remarkable degree 
of commonality in short-timescale kinetic properties. The distributions of transition-event 
durations tf, from all systems studied substantially "collapse" onto a single universal function 
after a simple linear rescaling (Fig. |J). Similarly, all the profiles of average local speeds 
during traversal of the reaction coordinate - - the "reaction-coordinate speed profiles" 
of the butane systems behave similarly (Fig. [5]). These profiles, furthermore, qualitatively 
match those from ID colored-noise systems, but are clearly distinct from the ID white-noise 
behavior (Fig. |6]). 

Despite this universal short-time behavior, quantitative estimates for reaction-rate, cor- 
relation, and event-duration timescales (Table p reveals no obvious relationship between 
butane systems with explicit solvents and those with implicit, stochastic "solvent." This 
suggests that stochastic simulations, at least at the level considered here, cannot be used 
for quantitative estimates of timescales in fully-solvated computations. While this assess- 
ment may be more pessimistic than that of Widmalm and Pastor ||, it is not surprising 
that consideration of an additional timescale, namely £&, makes agreement more difficult. 
Perhaps more important are two other issues: Do simple stochastic representations become 
better or worse as larger molecules are considered? What quantity of explicit solvent, say 
in combination with stochastic boundary conditions |27| , provides a good level of agreement 
with the array of kinetic quantities in fully-solvated systems? 

Further theoretical consideration of the connection between simple and molecular systems 
would also prove valuable. In particular, the study of simple models possessing a one- 
dimensional reaction coordinate coupled nontrivially to simple orthogonal degrees of freedom 
(based on, e.g., 0, ||, [l(| |ll|, |12| ) should permit a better understanding of the effects of noise 



color on the short-time transition behavior considered here. In a molecular system, as noted 



in Sec. |4.3| , the effective noise on any atom must be colored because of the coupling to the 
non-trivial "bath" of other particles. 

Following up on the motivation for the investigation reported here, we believe that an 
appreciation of the results can inform improved sampling approaches for reaction rate com- 
putations within the "transition path sampling" N22L |23| and "dynamic importance sampling" 



2TI O methods. These approaches require ensembles of transition-event trajectories. 



At a simple level, understanding the transition event duration £& in molecular systems is crit- 
ical: as noted in [Tj|, the duration of molecular crossing events sets, in part, the minimum 



computation time required to generate a realistic ensemble of transition-event trajectories. 
Furthermore, given the universal behavior uncovered here, stochastic simulations evidently 
reproduce key aspects of explicitly-solvated transition behavior. Despite the clear impor- 
tance of solvent effects, this study suggests that trajectories harvested in a computationally 
inexpensive stochastic context should — at a minimum — provide good starting points for 
refinement in a solvated context. 
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